#ifndef UHM_SOLVER_BASE_DIRECT_HXX_
#define UHM_SOLVER_BASE_DIRECT_HXX_

namespace uhm {

  class Direct_Solver_Base_ {
  public:
    Direct_Solver_Base_() { }
    virtual~Direct_Solver_Base_() { }

    virtual Int_ decompose()  { return false; }
    virtual Int_ updecompose()  { return false; }
    virtual Int_ solve()  { return false; }
    virtual Int_ check(Double_ &residual)   { return false; }
    virtual Int_ refine(Int_ niter, Double_ &residual,
                        Double_ thres=0.0) { return false; }
  };

  template<class Id_Nd_, class Id_El_, class Uhm_>
  class Direct_Solver_ : public Direct_Solver_Base_,
                         public Solver_<Id_Nd_,Id_El_,Uhm_> {
  public:

    // typedef is not inherited, explicitly bring Solver typeder here...
    typedef Solver_<Id_Nd_,Id_El_,Uhm_> _Solver_;

    typedef Mesh_<Id_Nd_,Id_El_,Uhm_> _Mesh_;
    typedef Tree_<typename _Mesh_::Element_*,UHM_TREE_MAX_CHILDREN> _Tree_;
    typedef Cache_<Id_Nd_,Uhm_> _Cache_;

    Direct_Solver_() 
      : Solver_<Id_Nd_,Id_El_,Uhm_>() { } 

    Direct_Solver_(Int_ type, Int_ uplo, Int_ sym,
                   Int_ rhs, Int_ blk)  
      : Solver_<Id_Nd_,Id_El_,Uhm_>(type, uplo, sym, rhs, blk) { }

    Direct_Solver_(Int_ type, Int_ rhs, Int_ blk) 
      : Solver_<Id_Nd_,Id_El_,Uhm_>(type, rhs, blk) { }

    //using Solver_<Id_Nd_,Id_El_,Uhm_>::get_flop;

    virtual~Direct_Solver_() { }

    virtual Int_ decompose() {
      time_in("uhm:: direct:: decompose",UHM_TIME_LEVEL_FRONT);        
#pragma omp parallel
      {
#pragma omp single nowait
        {
          for (Int_ i=0;i<this->tnds().size();++i) {
#pragma omp task firstprivate(i)
            {
              UHM_ERROR(op_leaf_to_me_par<_Solver_>(this, this->tnds()[i],
                                                    &op_uhm_direct_decompose<_Solver_>),
                        UHM_ERROR_DIRECT_DECOMPOSE);
            }
          }
        }
      }
      time_out(UHM_TIME_LEVEL_FRONT);

      return true;
    }

    virtual Int_ updecompose() {
      time_in("uhm:: direct:: updecompose", UHM_TIME_LEVEL_FRONT);
      this->cache().push_begin();

      time_in("uhm:: direct:: updecompose:: Up-Decomposing",UHM_TIME_LEVEL_INTERNAL);
#pragma omp parallel
      {
#pragma omp single nowait
        {
          for (Int_ i=0;i<this->tnds().size();++i) {
#pragma omp task firstprivate(i)
            {
              UHM_ERROR(op_leaf_to_me_par<_Solver_>(this, this->tnds()[i],
                                                    &op_uhm_up_direct_decompose<_Solver_>),
                        UHM_ERROR_DIRECT_DECOMPOSE);
            }
            
          }
        }
      }
      time_out(UHM_TIME_LEVEL_INTERNAL);

      time_in("uhm:: direct:: updecompose:: Upkeeping cache",UHM_TIME_LEVEL_INTERNAL);
      this->cache().push_end();
      time_out(UHM_TIME_LEVEL_INTERNAL);
      time_out(UHM_TIME_LEVEL_FRONT);
      return true;
    }

    virtual Int_ solve() {
      time_in("uhm:: direct:: solve", UHM_TIME_LEVEL_FRONT);
      time_in("uhm:: direct:: solve:: Phase 0",UHM_TIME_LEVEL_INTERNAL);

#pragma omp parallel
      {
#pragma omp single nowait
        {
          for (Int_ i=0;i<this->tnds().size();++i) {
#pragma omp task firstprivate(i)
            {
              UHM_ERROR(op_leaf_to_me_par<_Solver_>(this, this->tnds()[i],
                                                    &op_uhm_direct_solve_0_x<_Solver_>),
                        UHM_ERROR_DIRECT_SOLVE_ZERO);

              UHM_ERROR(op_leaf_to_me_par<_Solver_>(this, this->tnds()[i],
                                                    &op_uhm_direct_solve_1_x<_Solver_>),
                        UHM_ERROR_DIRECT_SOLVE_ONE);

              UHM_ERROR(op_me_to_leaf_par<_Solver_>(this, this->tnds()[i],
                                                    &op_uhm_direct_solve_2_x<_Solver_>),
                        UHM_ERROR_DIRECT_SOLVE_TWO);
            }
          }
        }
      }
      time_out(UHM_TIME_LEVEL_INTERNAL);
      time_out(UHM_TIME_LEVEL_FRONT);

      return true;
    }
    
    virtual Int_ check(Double_ &residual) {
      time_in("uhm:: direct:: check",UHM_TIME_LEVEL_FRONT);

      this->mesh().nodes2map(this->nodes());
      typename Uhm_::_B_ 
        b(this->get_type(), this->nodes().dof(), this->get_rhs()), 
        r(this->get_type(), this->nodes().dof(), this->get_rhs());

      b.create_flat_without_buffer(); b.create_flat_buffer();
      r.create_flat_without_buffer(); r.create_flat_buffer();

      time_in("uhm:: direct:: check:: Computing R = AX",UHM_TIME_LEVEL_INTERNAL);
      std::vector<typename _Cache_::Entry_*> w;
      this->work().dump(w);

#pragma omp parallel for schedule(static)
      for (Int_ i=0;i<w.size();++i) {
        Uhm_ uhm = w[i]->uhm();
        schur_x2r(uhm);

        Mapper_<Id_Nd_> s2nodes(&w[i]->schur(), &this->nodes());
        s2nodes.relate();
#pragma omp critical
        {
          axpy(s2nodes, uhm.BB(), b);
          axpy(s2nodes, uhm.RB(), r);
        }
      }
      time_out(UHM_TIME_LEVEL_INTERNAL);

      //std::cout << x << b;

      // b = b - r = b- Ax
      time_in("uhm:: direct:: check:: Computing B = B - R", UHM_TIME_LEVEL_INTERNAL);
      axpy(-1, r, b);
      norm1(b, residual);
      time_out(UHM_TIME_LEVEL_INTERNAL);

      b.free_all();                    
      r.free_all();

      time_out(UHM_TIME_LEVEL_FRONT);
      return true;
    }
    
    virtual Int_ refine(Int_ niter, Double_ &residual,
                        Double_ thres=0.0) {

      time_in("uhm:: direct:: refine", UHM_TIME_LEVEL_FRONT);
      
      this->mesh().nodes2map(this->nodes());
      typename Uhm_::_B_ 
        b(this->get_type(), this->nodes().dof(), this->get_rhs()), 
        r(this->get_type(), this->nodes().dof(), this->get_rhs());

      b.create_flat_without_buffer(); b.create_flat_buffer();
      r.create_flat_without_buffer(); r.create_flat_buffer();

      // Collects given matrices in work
      std::vector<typename _Cache_::Entry_*> w;
      this->work().dump(w);

      // Connectivity local2global
      time_in("uhm:: direct:: refine:: Connectivity local2global",UHM_TIME_LEVEL_INTERNAL);
      std::vector<Mapper_<Id_Nd_> > s2nodes; 
      s2nodes.reserve(w.size());
		
#pragma omp parallel for schedule(static)
      for (Int_ i=0;i<w.size();++i) {
        s2nodes.push_back(Mapper_<Id_Nd_>());
        s2nodes.back().set(&w[i]->schur(), &this->nodes());
        s2nodes.back().relate();
      }
      time_out(UHM_TIME_LEVEL_INTERNAL);

      // Iterative refinements
      for (Int_ iter=0;iter<niter;++iter) {
        time_in("uhm:: direct:: refine:: Computing Residual",UHM_TIME_LEVEL_INTERNAL);
        r.set_zero();
        b.set_zero();

        // R = AX
#pragma omp parallel for schedule(static)
        for (Int_ i=0;i<w.size();++i) {
          Uhm_ &uhm = w[i]->uhm();
          schur_x2r(uhm);

#pragma omp critical 
          {          
            axpy(s2nodes[i], uhm.BB(), b);
            axpy(s2nodes[i], uhm.RB(), r);
          }
        }

        // b = b - r = b- Ax
        axpy(-1, r, b);
        norm1(b, residual);
        time_out(UHM_TIME_LEVEL_INTERNAL);

        if (residual < thres)
          break;

        time_in("uhm:: direct:: refine:: Copying residual to UHM",UHM_TIME_LEVEL_INTERNAL);
#pragma omp parallel for schedule(static)
        for (Int_ i=0;i<w.size();++i) {
          copy_reverse(s2nodes[i], w[i]->uhm().RB(), b);
        }
        time_out(UHM_TIME_LEVEL_INTERNAL);

        if (get_verbose())
          std::cout << iter << " Iteration,  Residual = " << residual << "\n";

        // Residual is back to uhm but they are assembled values...
        // solve 0 needs to copy those values not merging
        time_in("uhm:: direct:: refine:: Solving for Error",UHM_TIME_LEVEL_INTERNAL);
#pragma omp parallel
        {
#pragma omp single nowait
          {
            for (Int_ i=0;i<this->tnds().size();++i) {
#pragma omp task firstprivate(i)
              {
                UHM_ERROR(op_leaf_to_me_par<_Solver_>(this, this->tnds()[i],
                                                      &op_uhm_direct_solve_0_r<_Solver_>),
                          UHM_ERROR_DIRECT_SOLVE_ZERO);
                UHM_ERROR(op_leaf_to_me_par<_Solver_>(this, this->tnds()[i],
                                                      &op_uhm_direct_solve_1_r<_Solver_>),
                          UHM_ERROR_DIRECT_SOLVE_ONE);
                UHM_ERROR(op_me_to_leaf_par<_Solver_>(this, this->tnds()[i],
                                                      &op_uhm_direct_solve_2_r<_Solver_>),
                          UHM_ERROR_DIRECT_SOLVE_TWO);
              }
            }
          }
        }
        time_out(UHM_TIME_LEVEL_INTERNAL);

        time_in("uhm:: direct:: refine:: Updating Solutioin",UHM_TIME_LEVEL_INTERNAL);
        for (Int_ i=0;i<w.size();++i) 
          axpy(1, w[i]->uhm().RB(), w[i]->uhm().XB());
        time_out(UHM_TIME_LEVEL_INTERNAL);
      }

      b.free_all();                    
      r.free_all();

      time_out(UHM_TIME_LEVEL_FRONT);
      return true;
    }
  };
}


#endif
